Stress field around arbitrarily shaped cracks in two-dimensional elastic materials 
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The calculation of the stress field around an arbitrarily shaped crack in an infinite two-dimensional 
elastic medium is a mathematically daunting problem. With the exception of few exactly soluble 
crack shapes the available results are based on either perturbative approaches or on combinations 
of analytic and numerical techniques. We present here a general solution of this problem for any 
arbitrary crack. Along the way we develop a method to compute the conformal map from the 
exterior of a circle to the exterior of a line of arbitrary shape, offering it as a superior alternative 
to the classical Schwartz-Cristoffel transformation. Our calculation results in an accurate estimate 
of the full stress field and in particular of the stress intensity factors Kj and Ku and the T-stress 
which are essential in the theory of fracture. 
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I. INTRODUCTION 

The existence of a crack in a stressed elastic medium 
does not necessarily mean that a catastrophic failure is 
in sight; otherwise many home owners would be in a con- 
stant state of panic. Indeed, one of the major objectives 
of the theory of fracture mechanics is to predict when the 
existence of a crack would necessarily lead to the failure 
of materials. A crucial ingredient of such a prediction is 
the determination of the state of deformation of a given 
material in the presence of the said crack under the ef- 
fect of a given external loading. This calculation is not 
available for arbitrarily shaped cracks even in two dimen- 
sional elastic media. To explain the difficulty consider for 
a moment an existing crack in an infinite two dimensional 
medium. Under given load conditions the displacement 
field u{r, t) is giving rise to the elastic strain tensor ey-: 



duj 



dxi 



(1) 



In linear elasticity the stress tensor is related to the strain 
tensor by 



E 



l-2v 



(2) 



where E and v are the Young's modulus and the Pois- 
son's ratio respectively. When the boundary conditions 
at infinity include both opening and shear modes (with 
respect to the straight crack) the stress field near the tip 
of any crack has a universal form i.e. 



Kr 



=4m 



(3) 



Here Kj and Kjj are the "stress intensity factors" with 
respect to the opening and shear modes, whereas J^jj{cp) 
and TifJ (ip) are universal angular functions common to 
all configurations and loading conditions. Despite the 
simplicity of Eq. JSJ the calculation of the stress intensity 
factors is not simple. Their numerical value depends on 



the shape of the crack, its dynamical history and on the 
far field boundary conditions. 

For a straight crack the stress intensity factors are 
known exactly j^'^l, i.e. 



Kr 



yy 
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(4) 



where 2a is the crack length and (t°° is the uniform load 
at infinity. The criterion for failure is then the famous 
Griffith- Irwin criterion [^, 



E 



(5) 



where F is the fracture energy. The physical meaning of 
Eq. is that the crack will initiate (and may cause 
failure) when the elastic energy flowing from the stress 
field in the bulk to the tip region is at least as large as 
the energy lost by lengthening the crack (bond breaking 
or any other energy cost involved). In the case of the 
straight crack failure will occur, for a given level of the 
external load, at a critical crack length. 

For cracks of arbitrary shape we still expect both Eqs. 
(O and (O to remain valid. The problem is then how 
to compute the stress intensity factors Ki and Ku when 
the crack is not straight. Indeed, a well known paper by 
Cotterell and Rice 6] offers such an analytic calculation 
for a slightly curved crack in perturbation theory in the 
amount of deviation from the straight crack. Many other 
works developed alternative hybrid numerical-analytical 
techniques like singular integral equations based either on 
dislocation distribution or crack "opening displacement 
functions" and superposition methods Q- A large body 
of research is devoted to direct numerical techniques like 
the finite element method 0- In this paper we offer a 
non-perturbative approach to the calculation of the full 
stress field of an arbitrarily shaped crack based on a con- 
formal mapping technique. In particular, we will be able 
to estimate accurately the stress intensity factors and any 
other relevant quantities like the T-stress, a quantity to 
be defined below (cf. Eq. 1401) '). In Sect. 2 we lay out the 
mathematical problem. In Sect. 3 we present a solution 
based on the method of iterated conformal maps. This 
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section includes also the general problem of finding the 
conformal map from the unit circle to a line of an arbi- 
trary shape. The section culminates in the calculation 
of the stress intensity factors. In Sect. 4 we exemplify 
the method and compare it against exactly soluble cases. 
Sect. 5 offers a short summary and a discussion. 



II. MATHEMATICAL FORMULATION 

The theory of elastostatic fracture mechanics in brittle 
continuous media is based on the equilibrium equations 
for an isotropic elastic body 



= 0. 



(6) 



For in-plane modes of fractures, i.e. under plane-stress 
or plane-strain conditions, one introduces the Airy stress 
potential U{x,y) such that 



dy 



2 ' '^^V 



dxdy 



dx"^ 



(7) 



Thus the set of Eq. lO, after simple manipulations, 
translate to a Bi-Laplace equation for the Airy stress po- 
tential U{x,y) 13 



AA[/(x,2;) = 0, 



(8) 



with the prescribed boundary conditions on the crack 
and on the external boundaries of the material. At this 
point we choose to focus on the case of uniform remote 
loadings and traction-free crack boundaries. This choice, 
although not the most general, is of great interest and 
will serve to elucidate our method. Other solutions may 
be obtained by superposition. Thus, the boundary con- 
ditions at infinity, for the two in-plane symmetry modes 
of fracture, are presented as 

axx{co) = ]ayy{co) = CToo ](yxy{oo) = Mode I (9) 
Oxxioo) = ;CTyy((X)) = ]axy{oG) = (Too Mode II . 

(10) 

In Addition, the free boundary conditions on the crack 
are expressed as 



O'xn(s) = 0-y„(s) = 



(11) 



where s is the arc-length parametrization of the crack 
boundary and the subscript n denotes the out-ward nor- 
mal direction. 

The solution of the Bi-Laplace equation can be written 
in terms of two analytic functions (j^iz) and ri{z) as 



[/(x,2/)=5RMz) + 77(z)] 



(12) 



In terms of these two analytic functions, using Eq. jT)), 
the stress components are given by 

ayy{x,y) = 3fi[2^'(z) + V(z)+?7"(z)] 
<7xx{x,y) = 3fi[2.^'(z) - z^"(z) - ?7"(z)] 
cjxy{x,y) = ^W{z)+ri"{z)]. (13) 



In order to compute the full stress field one should first 
formulate the boundary conditions in terms of the ana- 
lytic functions Lp{z) and 7]{z) and to remove the gauge 
freedom in Eq. \i'2.l . The boundary conditions Eq. Hll|) . 
using Eq. ((TJ, can be rewritten as ^ 



dt 



dU_ .dU_ 

dx dy 



= 



(14) 



Note that we do not have enough boundary conditions to 
determine U{x,y) uniquely. In fact we can allow in Eq. 
P2(l arbitrary transformations of the form 



ip ^ ip + iCz + 7 
V'-^'0 + 7, ip = rj' 



(15) 



where C is a real constant and 7 and 7 are complex 
constants. This provides five degrees of freedom in the 
definition of the Airy potential. Two of these freedoms 
are removed by choosing the gauge in Eq. I|14|l according 
to 



dx dy 



, on the boundary 



(16) 



It is important to stress that whatever the choice of the 
five freedoms, the stress tensor is unaffected; see 0] for an 
exhaustive discussion of this point. Computing Eq. I|lt)|) 
in terms of Eq. (|12|) we arrive at the boundary condition 



ip{z) + zip'{z) + tp{z) = 



(17) 



To proceed we represent ip{z) and tpiz) in Laurent ex- 
pansion form: 

ip{z) = ipiZ + ipQ + (p^l/z + ip^2/z'^ -\ , 

ij{z) = Viz-f +1/^-1/2 -fi/;_2/2^ + --- ■ (18) 

This form is in agreement with the boundary conditions 
at infinity that disallow higher order terms in z. One 
freedom is now used to choose ipi to be real and two 
more freedoms will allow us later on to fix ipo- Then, 
using the boundary conditions (j^ and (|10|l . we find 



"Pi 
^1=0 



^; ^1 = ^ Model, 



-01 = ido 



Mode II 



(19) 



III. THE SOLUTION 



As said above, the direct determination of the stress 
tensor for a given arbitrary shaped crack is difficult. To 
overcome the difficulty we perform an intermediate step 
of determining the conformal map from the exterior of 
the unit circle to the exterior of our given crack. Cur- 
rently the best available approach for such a task is the 
Schwartz-Cristoffel transformation. Here we will present 
an alternative new approach for finding the wanted con- 
formal transformation, given in terms of a functional it- 
eration of fundamental conformal maps. The use of it- 
erated conformal maps was pioneered by Hastings and 
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The inverse mapping 



'0=0, A 



is of the form 






FIG. 1; Example of how to construct the conformal mapping 
along a line. 



Levitov it was subsequently turned into a powerful 
tool for the studyof fractal and fracture growth patterns 
ES in m El 111 111 El- In the next subsection we 
describe how, given a crack shape, to construct a con- 
formal map from the complex cj-plane to the physical 
z-plane such that the conformal map z = <i>(w) maps 
the exterior of the unit circle in the oj-plane to the exte- 
rior of the crack in the physical z-plane, after n directed 
growth steps. We draw the reader's attention to the fact 
that this method is more general than its application in 
this paper, and in fact we offer it as a superior method 
to the Schwartz-Cristoffel transformation, with hitherto 
undetermined potential applications in a variety of two- 
dimensional contexts. 



A. The conformal mapping 

The essential building block in the present application, 
as in all the applications of the method of iterated con- 
formal maps is the fundamental map that maps the 
exterior circle onto the unit circle with a semi-circular 
bump of linear size \f\ which is centered at the point 



e . This map reads [Qj: 



2w 



1 + w + w [ 1 + 



2 1- A 

wl + X 



1/2 



(20) 

1/2 

(21) 



\Z ~ VTTA(z2 - 1) 

1 - (1 + A)z2 ' 



(22) 



By composing this map with itself n times with a judi- 
cious choice of series {9k}k=i and {Xk}n=i "^'^^ con- 
struct $("^(a;) that will map the exterior of the circle 
to the exterior of an arbitrary simply connected shape. 
To understand how to choose the two series {Ok}^^i and 
{Xk}n=i consider Fig. ^ and define the inverse map uo — 
x'^")(z). Assume now that we already have <I>*^"~^^ (w) and 
therefore also its analytic inverse X^"~^''(^) after n — 1 
growth steps, and we want to perform the next iteration. 
To construct $'■"'(0;) we advance our mapping in the di- 
rection of a point z in the z-plane by adding a bump in 
the direction oi w — X^"^^''(-S) in the w-plane. The map 
<!)(") (w) is obtained as follows: 



$(»)(C^) =$("-!) (^,^^ _;, J,,)) 

The value of 0„ is determined by 

0„ = arg[x("-^H5)] 



(23) 



(24) 



The magnitude of the bump A„ is determined by requir- 
ing fixed size bumps in the z-plane. This means that 



Xn — 



An 



(e' 



(25) 



We note here that it is not necessary in principle to have 
fixed size bumps in the physical domain. In fact, adap- 
tive size bumps could lead to improvements in the pre- 
cision and performance of our scheme. We consider here 
the fixed size scheme for the sake of simplicity, and we 
will show that the accuracy obtained is sufficient for our 
purposes. Iterating the scheme described above we end 
up with a conformal map that is written in terms of an 
iteration over the fundamental maps (|20(l : 



$(")(W) = O . . . O (/)e„a„(w) • 



(26) 



For the sake of newcomers to the art of iterated conformal 
maps we stress that this iterative structure is abnormal, 
in the sense that the order of iterates in inverted with 
respect to standard dynamical systems. On the other 
hand the inverse mapping follows a standard iterative 
scheme 



(27) 



The algorithm is then described as follows; first we 
divide the curve into segments separated by points {zi\. 
The spatial extent of each segment is taken to be approx- 
imately VAo , in order to match the size of the bumps in 
the z-plane. Without loss of generality we can take one 
of these points to be at the center of coordinates and to 
be our starting point. From the starting point we now 
advance along the shape by mapping the next point Zi 
on the curve according to the scheme described above. 



4 



B. Solution in terms of conformal mappings 

The conformal map $'^"^(0;) is constructed in n itera- 
tive steps. For the discussion below we do not need the n 
superscript and will denote simply ^{uj) = This 
map is univalent , having the Laurent expansion form 



(28) 



Any position z in the physical domain is mapped by 
x{z) = $^^(z) onto a position uj in the mathematical 
domain. This transformation does not immediately pro- 
vide the solution as the Bi-Laplacian operator, in con- 
trast to the Laplacian operator, is not conformally in- 
variant. Nevertheless, the conformal mapping method 
can be extended to non-Laplacian problems. We begin 
by writing our unknown functions (p{z) and '>jj{z) in terms 
of the conformal map 



viz) = 'P ix{z)) , ^{z) = -ip ixiz)) 



(29) 



Using the Laurent form of the conformal map, Eq. I|28|l . 
the linear term as a; ^ cx) is determined by Eqs. (|29|1 . 
We therefore write 

ip{uj) = ipiFiu; + ip^i/u; + (f^2/^'^ + ■ ■ ■ , 

Vi(w) = Vi-Fit^ + -0o + V'-i/^ + V'-2/t^^ + --. ,(30) 

where we used the last two freedoms to choose ipo = 
—Fqcpq such that ipo — 0. The boundary condition 117|) 
is now read for the unit circle in the cj-planc. Denoting 



u{e) 



00 

E 

Tl = l 



ip-^n/e" , v{e) 



n=0 



we write 



uie) 



$(e) 

$ (e) 



(31) 



(32) 



The function /(e) is a known function that contains all 
the coefficients that were determined so far: 



we need to Fourier transform the function $(e)/$'(e). 
This is the most expensive step in our solution. One 
needs to carefully evaluate the Fourier integrals between 
the branch cuts. Using the last two equations together 
with Eqs. (EU and ^ we obtain 

00 

fc=l 
00 

V"*™ - '^k bm-k-i'flk = fm ■ m 0, 1,2 • • • (37) 

k=l 

These sets of linear equations are well posed. The coef- 
ficients ^-m can be calculated from equation (|36|l alone, 
and then they can be used to determine the coefficients 
i/j-m. This is in fact a proof that Eq. H32(l determines 
the functions u(e) and v{e) together. This fact had been 
proven with some generality in 

The calculation of the Laurent expansion form of i^(w) 
and ipi^) provides the solution of the problem in the co- 
plane. Still, one should express the derivatives of ip{z) 
and ri{z) in terms of <^(w) and iP{lo) and the inverse map 
xiz) to obtain the solution in the physical z-plane. This 
is straight forward and yields 

^'(z) = ^'[x{z)] x'{z) 

^"{z) = ^"[x{z)] [X\z)? + V'[x{z)]x"{z) 
r^"{z) = ^'(z) = ^'[x(z)] ^{z). (38) 

Upon substituting these relations into Eq. Ijl^fl one can 
calculate the full stress field for an arbitrarily shaped 
crack. The expression of the stress field in terms of the 
inverse conformal mapping is known for quite a long time 
although it is very limited as the conformal mapping and 
its inverse is rarely at hand. The central step of progress 
in this paper is the conjunction of the novel functional 
iterative scheme for obtaining the inverse conformal map- 
ping with the known result that expresses the stress field 
in terms of this inverse mapping. 



/(e = -ipiFie ~ ==ipiFi 

$'6 e 



(33) 



To solve the problem we need to compute the coefficients 
(fin and ^fjn- To this aim we first write |T5l | 



$'(e) 



(34) 



The function /(e) has also an expansion of the form 



(35) 



C. The Stress Intensity Factors and the T-stress 

Linear elasticity fracture mechanics, under small scale 
yielding, has no intrinsic length scale. Accordingly the 
stress intensity factors are the most important quanti- 
ties that characterize the universal near-tip fields. At 
this point we explain how to calculate the stress inten- 
sity factors from our solution. 

In principle, the calculation follows directly from the 
solution in terms of the conformal map as described 
above. Previous authors derived the following expression 
for the complex combination (of the real) stress intensity 
factors 



In the discussion below we assume that the coefficients bi 
and fi are known. In order to compute these coefficients 



Ki - iKii = 2 



Ifi'iujj), 



(39) 
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where Uj is the position of the tip in the w-plane, and 
5j is the argument of the tip position Zj in the physi- 
cal plane. This result, although exact, cannot be used 
to obtain accurate estimates of the stress intensity fac- 
tors in our method. Since we construct our crack from a 
succession of little bumps, our crack tip is not infinitely 
sharp. Therefore we cannot base our calculation of the 
stress intensity factors on the precise coordinates of the 
tip. Rather, we need to exploit our knowledge of the 
stress field for a substantial region around the tip. We 
thus consider the stress field in the region of the tip, and 
write the components along the tangent to the crack at 
the tip 




(40) 



We have added terms of 0{^) to the leading terms, and 
in addition we took into account the T-stress contribu- 
tion to the purely radial component (!„■ It is crucial to 
note that there are no other 0(1) terms in the first two 
lines of Eq. H40(l . We did not need to consider explicitly 
any higher 0(r'^/^) terms. In order to extract the stress 
intensity factors from the full field distributions we fit 
our solution near the tip to the form given by Eq. (|40|l . 
We thus obtain not only the stress intensity factors but 
also the coefficient of the subleading terms. We note that 
the so-called T-stress has an important role in fracture 
theory (cf. 0,1111). 




FIG. 2: The stress tensor component Gyy for a straight crack 
with a = 200. The component is evaluated for y — 
and a distance p away from the tip at a; = 200. The nu- 
merical estimates 5-yy approaches the analytical result (the 
left most line) as we decrease the linear size of the bumps 
\/Ao = 1.0, 0.9, 0.8, . . . , 0.3. The inset shows the relative er- 



^f^TT^ (z2 - a2)3/2 



(^2 _ a2)3/2 



(41) 



IV. DEMONSTRATIONS OF THE METHOD 

In this section we demonstrate our method for two 
(of the very rare) exactly soluble geometries. There are 
many numerically solved problems in the literature that 
one can use for comparison, but our aim here is to ex- 
emplify the essentials of our approach. We will see that 
the method works very well and we propose that it can 
be used for arbitrary crack shapes as well. 



A. The Straight Crack 

The problem of a crack of length 2a in an infinite do- 
main, subjected to a remote uniaxial load Gyy = a°° and 
traction-free crack faces is considered as the canonical 
problem in the theory of linear elasticity fracture me- 
chanics. The resulting stress field is known, and is given 
by i: 



(Tyy{x,y) = cr°°5R 



tya^ 



(z2 



2)3/2 



where z =^ x + iy and the crack is represented by a branch 
cut along —a < a; < a, y — 0- 

We applied our general solution to the straight crack 
problem. We constructed the conformal mapping using 
our functional iteration scheme although the exact con- 
formal mapping is known to be 



<i>(.) = |L+l). 



(42) 



Fig. 121 compares our calculation of ayy along the x-axis 
with the exact result. The deviation near the tip of the 
crack is expected as in our solution the crack tip is not 
infinitely sharp but has a finite radius of curvature con- 
trolled by the parameter A in Eq. 1)20(1 . Decreasing the 
value of A, we can obtain more accurate results. Note 
that real crack-tips do have a finite radius of curvature, so 
that our method can be even more appropriate; the ide- 
alization of representing cracks by mathematical branch 
cuts is no longer a necessary simplification. 
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FIG. 3: The stress tensor components along the tangent to 
the crack tip. The crack is a semi circular arc of radius 200. 
p is the distance away from the tip. The points show the 
numerical values computed using bumps of linear size \/Ao = 
0.6 and the lines are the corresponding analytical values. 



FIG. 4: Stress intensity factors for circular arcs of angles 29. 
The points are the numerical values computed from fitting of 
the stress field along the tangent to the crack tip. The fitting 
function used is of the form Eq. 1401 . In this calculation 
y/Xo = 0.4, and the fitting window is 5 < p < 30. 



B. The Circular Arc Crack 

The analytic methods developed by Muskhelishvili 
lead to the solution of various problems concerned with 
circular regions and infinite regions cut along circular 
arcs. Here we consider a crack in the shape of a circular 
arc that extends from z = e"^^ to z — e~'^. This crack 
is subjected to a remote uniaxial load a^x — parallel 
to the X-axis with traction-free crack boundaries. The 
stress tensor components are calculated from Eq. (|13f) 

with Hi 



-t- 

n{z) = 

^'(z) 
Here 



2 G(z) 



C{d){z - cos{9)) 



1 cje) 

4 2 



1 

4^2 



2 G{z 

a. 



C{9){z ~ cos{0)) 



cos{e) 
2z 



cos{e) 

2z 



1 

2^2 



1 

2^2 



1 Cje) 

4 2 



1 

4^ 



(43) 



c{e) = 

G{z) = 



1 - sin^ (6'/2)cos2 [9/2) 

2[1 + sin^ {0/2)] 
v/z2 - 2zcos(0) + 1, 



(44) 



where G{z) is defined such that z^^G{z) ^ 1 as z ^ oo, 
which leads to G(0) = — 1. Note that these results are 
given for an arc of unit radius and therefore the solution 
for any other circular arc can be obtained by a suitable 
rescaling transformation. 

We applied our method for this configuration. The 
results for all the three components of the stress tensor 
field for = 7r/2 along the continuation of the tip parallel 
to the X-axis are presented in Fig. |2| and compared with 
the exact results. 

Finally, we used the method to extract the stress in- 
tensity factors. Fig. 0] shows the mode I and mode II 
stress intensity factors as a function of the half arc an- 
gle 9 and compares them with the exact analytic result 
derived form the full field solution 



V. SUMMARY AND CONCLUSIONS 

In summary, we have demonstrated that the method 
of iterated conformal maps can be used to construct the 
conformal map from the exterior of the unit circle to the 
exterior of an arbitrary crack. Next one can address the 
calculation of the stress field around such a crack with 
given loads at infinity. Having the said conformal map 
at hand simplifies enormously the calculation of the full 
stress field, allowing an accurate estimate of the stress 
intensity factors or of the sub-leading terms like the T- 
stress. The method was demonstrated by comparison 
with exactly soluble examples. The quality of the com- 
parison leads to the conclusion that the issue of potential 
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failure of a material given a crack and boundary condi- 
tions can be efficiently dealt with. Future work will em- 
ploy the present method to describe the dynamics of slow 
fracture where quasi-static methods are adequate. 
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